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Abstract 

The thermodynamic, dynamic and structural behavior of a water-like system confined in a ma¬ 
trix is analyzed for increasing confining geometries. The liquid is modeled by a two dimensional 
associating lattice gas model that exhibits density and diffusion anomalies, in similarity to the 
anomalies present in liquid water. The matrix is a triangular lattice in which fixed obstacles im¬ 
pose restrictions to the occupation of the particles. We show that obstacules shortens all lines, 
including the phase coexistence, the critical and the anomalous lines. The inclusion of a very dense 
matrix not only suppress the anomalies but also the liquid-liquid critical point. 
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I. INTRODUCTION 


The phase behavior of systems as particles interacting via the so-called core-softened 
(CS) potentials has received a lot of attention recently. They show a repulsive core with a 
softening region when particles are very close and an attractive region when particles are 
more distant. These CS can be modeled as continuous potentials or lattice gas models. For 
the lattice structure the two competing scales arise from two equilibrium conhgurations: 
low density and high density. This procedure generates models that are analytically and 
computationally tractable and that one hopes are capable of retaining the qualitative features 
of the real complex systems. The physical motivation behind these studies is the assumption 
that two length scales systems exhibit the same anomalous behaviors present in water. 

Confirming this hypothesis a number of continuous [1, 2, 4-7] and lattice gas models [8- 
16] show the presence of density, diffusion and structural anomalous behavior as observed 
in water [17, 18]. 

In addition to the seventy-two anomalies [19], water has at very low temperatures two 
coexisting amorphous phases with distinct densities: the low density amorphous (LDA) and 
high density amorphous phases (HDA) [20-24]. These two amorphous phases led to the 
hypothesis of the existence at higher temperatures of two liquid phases: a low density liquid 
and high density liquid phases. Such conjecture establishes that the coexistence between 
these two liquid phases ends in a second critical point or also called, liquid-liquid critical 
point (LLCP) [24]. Experiments for testing the existence of this criticality are difficult 
since the region in the pressure versus temperature phase diagram where the alleged criti¬ 
cal point exists is locate beyond the homogeneous nucleation limit. In order to circumvent 
this difficulty for testing the existence of the liquid-liquid critical point recently confined 
geometries have been employed [25, 26]. In these nanoconfined geometries the disruption 
of the hydrogen bonds suppress the solidification of the system and allows for maintaining 
the system liquid in temperatures in which otherwise would be solid [25, 26]. These experi¬ 
mental systems show convincing evidences that water exhibits two liquid structures at low 
temperatures. 

The use of confining water, however brings another set of issues. What guarantees that 
the same thermodynamic and dynamic anomalies and criticality present in the bulk are 
not destroyed as the system is confined? Can confinement bring up new phenomena not 
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present in the bnlk system? In order to answer these qnestions water-like atomistic or 
continuous effective potential models were explored confined by different geometries such as 
plates [27-33], one pore [34-42] and disordered matrix [39, 48-52], These simulations show 
that confinement leads to a controversial result regarding the melting line. While results for 
SPC/E water show that the melting temperature for hydrophobic plates is lower than the 
bulk and higher than for hydrophilic walls, for mW model no difference between the melting 
temperature due to the hydrophobicity [43] is found. 

In confined systems the TMD occurs at lower temperatures for hydrophobic confine¬ 
ment [44, 45] and at higher temperatures for hydrophilic confinement [46] when compared 
with the bulk. The diffusion coefficient, D, in the direction parallel to the plates exhibit an 
anomalous behavior as observed in bulk water. However the temperatures of the maximum 
and minimum of D are lower than in bulk water [45]. In the direction perpendicular to the 
plates, no diffusion anomalous behavior is observed [47]. 

In addition to the usual density and diffusion anomalous behavior, these confined systems 
show a variety of new effects not present in the bulk. For example, fluids confined in carbon 
nanotube exhibit formation of layers, crystallization of the contact layer [53, 54] and a 
superffow not present in macroscopic confinement [55, 56]. 

The confinement by a pore, within plates or nanotubes is symmetric and even though 
it introduces a breaking of the water hydrogen bond network, this is done in an ordered 
way. Confining matrix such as the ones present in plants and underground water are not 
ordered. Recently the effects of an water-like liquid confined in a disordered matrix have 
been analyzed using a model in which the bonds are introduced by the inclusion of Potts 
variables [49]. This study shows that the liquid-liquid coexistence line is affected by the 
increase of the density of random porous in the matrix without disappearing. 

In all these studies, however the liquid-liquid transition is preserved and the density of 
confining matrix is not very high. Here, we give a further step by investigating effects 
imposed by disordered porous when the random matrix exhibits a high density. The system 
is defined in a triangular lattice where the obstacles are fixed and randomly distributed. The 
fluid is modeled as a Associating Lattice Gas Model [12] defined an occupational variable 
together with a bond orientational variable. This model in the bulk shows the density 
and diffusion anomalies present in water and the liquid-gas and liquid-liquid criticality [12], 
Here we explore the effect in the chemical potential versus temperature phase diagram of 
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the presence of the random hxed obstacles. 

This paper is organized as following. In Sec. II we present the model used here. In Sec. 
Ill outline details of Monte Carlo simulations and how was calculated the thermodynamic 
properties of system. Results are presented in Sec. IV, followed by conclusions in Sec. V. 


II. THE MODEL 


The model system is dehned in a triangular lattice. Each accessible site i can be empty 
or occupied by a water molecule. Empty sites have Uj = 0 while occupied sites have (Tj = 1. 
Each water molecule has orientational states represented by the variable r that presents six 
arms, being two inert arms with r* = 0 and four active arms with Tj = 1. They represent 
the possibility of a molecule to form hydrogen bonds with up to four neighbor molecules. 
For representing the symmetry present in water, two inert arms are diametrically positioned 
and just three different orientational states are possible. Fig. 1 exemplihes the geometry of 
a water molecule. A hydrogen bond is formed only when the active arms of two neighbor 



FIG. 1. The occupational and orientational states of a water molecule placed at the site i. In such 
example, the arms variables of molecule read: = 0, = 1, = 1, = 0, = 1, 

rr = 1 . 

molecules point out to each other, TiTj = 1. In this case, the interaction energy between 
two bonded arms reads —v and while non bonded arms contribute with a higher energy of 
—V + 2u (punishement for non forming hydrogen bonds). The Hamiltonian of the system is 
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given by 


( 1 ) 


Ti = 2u (TjCTj 

(*.i> 




The phase behavior of the system in the absence of obstacles was already analyzed [15]. 
Below it will be reviewed. At ground state, T* = T/v = 0, the grand potential per site 
is $ = e — fiN where e = (H) /T^. For low chemical potentials, the lattice is empty 
and the system is constrained in gas phase, p = 0. In this phase the grand potential 
is ‘Fgas = 0. Increasing the chemical potential the system reaches a point at which the 
gas phase coexists with a low density liquid phase (LDL). In this phase, the density is 
p = 3/4: and each particle forms four hydrogen bonds with its neighbor, resulting in a grand 
potential per site ^ldl/L"^ = ~(3/2)n — (3/4)/r and consequently in a gas-LDL coexistence 
chemical potential Pg-ldl = f^G-LDh/v = —2. For high chemical potentials, all sites of 
lattice are occupied by particles, resulting in a density p = I and grand potential per site 
^hdl/L'^ = —3n + 2u — p. The coexistence between the LDL phase and the HDL phase 
occurs at Pldl-hdl — I^ldl-hdl/v = Su/v — 6. The main features of LDL and HDL 
phases are exemplihed in Fig. 2 for two possible conhgurations at T* = 0. At temperatures 



FIG. 2. (a) HDL phase for the bulk system. The solid lines indicate the bonding arms (b) LDL 
phase for the bulk system. The solid lines indicate the bonding arms. 

T* > 0, the model was studied by Monte Carlo simulations whose phase diagram is shown 
in Fig. 3. The gas-LDL and LDL-HDL transition lines are hrst-order transitions, ending in 
respective tricritical points Td and T^, respectively that are joined by a line of continuous 
transitions, called A—line. For the bulk case, T^i and read 0.65 and 0.825, respectively. 
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FIG. 3. For the bulk System, panel (a) shows the phase diagram /r* vs T*, illustrating the gas- 
LDL (empty circles) and the LDL-HDL (filled circles) phase transitions, the A-line (empty squares) 
and the TMD line (filled triangles). In panel (b) we plot the Cy versus T* for = —0.80 (circles), 
fj,* = 0.60 (diamonds) and fi* = 1.20 (squares). In (c) the system density p versus T* for fixed 
along the TMD line (dashed line). 


In order to understand the differences between the LDL and HDL phases, we divide the 
lattice in four sublattices as illustrated in the Fig. 4. The LDL phase is characterized by 
one of the sublattices being empty while all the others are hlled, in such a way that the 
transition to the HDL phase occurs when the empty sublattice is hlled. Also, it is signed by 
a rotation in the inert arms, in which in the HDL phase they are all parallel. In addition, in 
the LDL phase each particle forms four bonded arms that show a zigzag structure, whereas 
in the HDL phase each particle also forms four bonded arms but in addition to the zigzag 
structure two parallel lines appear. In the regime of very high temperatures, the system is 
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disordered, in which the sublattice occupations do not exhibit any ordering. By lowering T, 
the A-line is crossed, which one sublattice is emptied and the others remaining filled with 
an reorganization of the inert arms that form the above ordered zig-zag structure. 

The density of bonds, phb = yy J2i=i X^i +5 is also an important quantity for 

characterizing the phase transitions. At T* = 0 the gas, LDL and HDL phases has phb 
reading 0,1.5 and 2, respectively. Thus the phase transitions are also signed by changes 
in the fraction of hydrogen bonds. At high temperature the system is disordered and the 
A—line is obtained through the specific heat at constant volume cy by 


Cy = 


1/T2 


\ / e 




+ 


3Nk 


B 


2V 


( 2 ) 


where 6X = X — (X) with X = Ti and N and averages are evaluated in the ensemble of 
T, p fixed. 

In this work the porous matrix is introduced by considering fixed obstacles that are 
randomly distributed in the lattice. Each obstacle occupies a single site and interacts with 
the particles via a “hard core” constraint. The density of obstacles is given by po = No/L‘^ 
where No is the number of obstacles and is the system volume. In Fig. 4, we exemplify 
a lattice configuration composed of water, obstacles and empty sites. 



FIG. 4. Example of a lattice configuration. Gray, black and dashed circles describes water 
molecules, obstacles and empty sites, respectively. For clarity, hydrogen bonds have not explicitly 
shown. 
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III. THE METHODS AND SIMULATION DETAILS 


All the thermodynamic properties have been obtained by performing grand canonical 
Monte Carlo (MC) simulations for hxed and po [57]. Microscopic conhgurations are 

generated according to the Metropolis algorithm [58] described as follows. Obstacles are 
initially randomly distributed and a given site k other than a porous is randomly chosen. 
Two sorts of transitions are possible, taking into account k be empty or occupied by a water 
molecule. In the former case, a water molecule in one of its arm states is chosen, whereas 
in the latter one of the three possibilities (including the other two arm states and empty 
site) are performed. All of possible transitions are chosen with equal probability. Next, we 
evaluate the energy difference ATf between the original and the new conhguration. The 
conhguration change is accepted according to the Metropolis prescription min{l, 
where (3 = l/fc^T. A Monte Carlo step is dehned as the number of trials for generating new 
conhgurations. After discarding a sufficient number of initial MC steps, relevant quantities 
are evaluated. In addition to the thermodynamic quantities, we also investigated the inhu- 
ence of obstacles in the dynamic properties, characterized through the diffusion coefficient 
D given by Einstein’s relation 


D 


lim 

t^OO 


(Ar(t)^) 

At 


( 3 ) 


where (Ar(f)^) = ((r(f) — ^(0))^) is the mean square displacement per particle and time is 
measured in Monte Carlo steps. The numerical MC procedure for calculating the diffusion is 
described as follows. First, the system is equilibrated by employing the previous Metropolis 
dynamics for hxed T* and p*. After the equilibrium is reached, an occupied site i and it’s 
neighbor j are chosen randomly. In case of neighbor site j be empty, the molecule moves to 
the empty site also following the above Metropolis prescription min{l, where ATi 

is the diherence of energy due to the movement. A Monte Carlo step is dehned through the 
number of trials of movement for every particle. After repeating this algorithm nt times, 
where n is the number of molecules in the lattice, the dihusion coefficient is calculated from 
Eq.(3). 

Numerical simulations have been performed for triangular lattices of size L = 35 and 
periodic boundary conditions for three representative values of density of obstacles po = 
0.08, 0.24 and 0.40 have been considered. In all cases, we have used 10® Monte Carlo (MC) 
steps to equilibrate the system and 10® MC steps for evaluating the relevant quantities. 



IV. RESULTS 


A. Structural and thermodynamic behavior 





FIG. 5. p vs p* for distinct porous densities po for T* = 0.40. 


First, let us exam what happens with the phases present in the bulk system as the 
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obstacles are introduced. Fig. 5 shows the water density p, versus the reduced chemical 
potential p*, for distinct porous densities at the hxed temperature T* = 0.40. The inclusion 
of obstacles changes the gas-LDL and the LDL-HDL phase transition, whose effect becomes 
more pronounced as po increases. In particular, by raising the porosity, the density gap 
between the liquid phases becomes less abrupt and inclusion of obstacles move the transition 
points for larger chemical potentials. Figs. 6, 7 and 8 illustrate the chemical potential 
versus temperature phase diagrams for po = 0.08, 0.24 and 0.40, respectively. In particular, 
by increasing p^ the tricritical points T^i and Tc 2 , in which the gas-LDL and LDL-HDL 
coexistence lines meet the A—line, decreases as shown in Figs. 6 , 7 and 8. More specihcally, 
while the bulk gas-LDL tricritical point is located at Td = 0.65, it moves to Td = 0.60, 0.55 
and 0.52 for po = 0.08, 0.24 and 0.40, respectively. 

This scenario becomes even more drastic in the case of the LDL-HDL phase transition. 
The tricritical point not only decreases its value from = 0.825 (bulk) to Tc 2 = 0.57 and 
Tc 2 = 0.52 for Po = 0.08 and 0.24, respectively but it disappears for po = 0.40, implying 
the absence of liquid-liquid transition line. In addition the A-line can also be found only 
for po = 0.08, 0.24, in which for po = 0.40 the original bulk tricritical Tqi becomes a single 
critical point . 

The above changes in liquid phases as well as transition points can be understood by 
verifying that the inclusion of porous suppress partially the structured patterns found in 
the LDL and HDL phases (see e.g Fig. 2(a) and (b) for the bulk case). In the case of the 
LDL phase the ordered bulk structure is distorted as po increases, as illustrated in the Fig. 
9 for p* = —0.5. For the lowest case po = 0.08, the degree of conhnement is low and most 
occupied sites preserve at least three bonds. As the density of obstacles is increased (here 
exemplihed for po = 0.16 and 0.24) the fraction of disrupted bonds increases, reaching a 
limit in which the connectivity of the network is lost. Similar effect is verihed in the HDL 
phase, but the effect is more pronounced in such case. This can be understood by recalling 
that in contrast to the HDL phase, porous occupy partially empty sites with neighboring 
molecules not forming hydrogen bonds in the LDL phase. This lost of connectivity also 
explains, in similarity with the decrease of tricritical points, why the transition from the 
disordered structure to the LDL through the A-line occurs for lower temperatures than the 
bulk. Recalling that such transition is characterized by the disordered phase ordering by 
making one of the sublattices empty, the inclusion of porous makes the system entropy 
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FIG. 6. For po = 0.08, in (a) the phase diagram in the plane of reduced chemical potential p* versus 
reduced temperature T*. Empty and filled circles denote the Gas-LDL and the LDL-HDL phase 
transitions, respectively. The A and TMD lines are described by empty squares filled triangles, 
respectively. Panel (6) shows the specific heat at constant volume cy versus T* for the system with 
obstacles (filled squares) and bulk system (empty circles) at /r = —1.00. Panel (c) shows p versus 
T* for p* = 0.0,..., 2.2 showing the TMD line (dashed line). 
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be larger than the bulk, with partial disruption of hydrogen bonds. Thus, all transition 
points, move for lower temperatures as a way for ’’compensating” the above increase of 
disorder. In other words, due to the inclusion of porous, the structured phases exist only 
for lower temperatures than the bulk, whose decreasing become more pronounced as po 
increases. Finally, for high density of obstacles the transition is destroyed by enhancement 
of fluctuations. The last comment concerns in the comparison between the TMD as po 
increases. As for the transition lines, the TMD shortens and move for lower temperatures 
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FIG. 7. For po = 0.24, in (a) the phase diagram in the plane of reduced chemical potential p* versus 
reduced temperature T*. Empty and filled circles denote the Gas-LDL and the LDL-HDL phase 
transitions, respectively. The A and TMD lines are described by empty squares filled triangles, 
respectively. Panel (6) shows the specific heat at constant volume cy versus T* for the system with 
obstacles (filled squares) and bulk system (empty circles) at /r = —1.00. Panel (c) shows p versus 
T* for p* = 0.0,..., 2.0 showing the TMD line (dashed line). 

(with maximum p decreasing) as po increases. However, in contrast previous results, for 
Po = 0.40 a tiny TMD (ranged from T* = 0.50 to 0.70 with p = 0.56 to 0.64) is verified. 


B. Diffusion and dynamic anomaly 

Besides the influence of immobile obstacles in the thermodynamic quantities, another 
relevant question concerns what happens with the system mobility as the density of obstacles 
increases. Fig. 11 shows the diffusion coefficient computed using Eq. (3) for different T*’s 
and po’s. In similarity with the pure model, the diffusion coefficient presents an increasing 
anomalous behavior until a maximum value by raising p for po = 0.08 and 0.24. The only 
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FIG. 8. For po = 0.40, in (a) the phase diagram in the plane of reduced chemical potential p* versus 
reduced temperature T*. Empty and filled circles denote the Gas-LDL and the LDL-HDL phase 
transitions, respectively. The A and TMD lines are described by empty squares filled triangles, 
respectively. Panel (6) shows the specific heat at constant volume cy versus T* for the system with 
obstacles (filled squares) and bulk system (empty circles) at /r = —1.00. Panel (c) shows p versus 
T* for p* = 0.2,..., 2.2 showing the TMD line (dashed line). 

exception is the case with density of obstacles 0.40. 

The reason concerns that the dynamic anomaly depends crucially of the presence of a 
high number of neighbor sites occupied by the fluid [62]. The obstacles make this difficult 
and for a very high number of obstacles, the mobility becomes even impossible. 

Since for water-like systems, typically the region in the p*-T* phase diagram in which 
the density anomaly is present is close to the region where the diffusion anomaly appears. 
Therefore one expects that the suppression of the hrst is directly related to the disappearance 
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FIG. 9. Spatial snapshot (35 x 35 sites) of triangular lattice. Each site is represented by hexagon, 
with its six nearest-neighbor sites. White hexagons represent vacancies, black represent obstacles 
and gray represent water-like particles. The snapshots exhibit character configurations of system 
with chemical potential fi* = —0.5 and temperature T* = 0.3. In (a) we present the bulk system. 
In (b) the system submitted at low degree of confinement po = 0.08 and the blue rectangles denote 
the regions where the characteristic geometry of LDL of ALG is preserved. In (c), intermediate 
degree of confinement p = 0.24, and green rectangles denote the LDL structure. The highest degree 
of confinement p = 0.40 is shown in (d). 

of the other. 


V. CONCLUSION 


The effects of hxed obstacles in thermodynamic and dynamic properties of an simplihed 
water-like model have been investigated. For low degree of conhnement, the thermodynamic, 
structural and dynamic properties of model are almost totally preserved due to the low steric 
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FIG. 10. Chemical potential versus temperature illustrating (a) the TMD lines (b) the gas-LDL 
tricritial point (c) the LDL-HDL critical point values for the bulk and the system with different 
concentrations of obstacles. 


effects. For intermediate case, po = 0.24, the system suffers significant changes such as, the 
decrease of the critical and tricritical points to lower temperatures, resulting in a reduction 
of coexistence regions. This effect is more dramatic for the liquid-liquid coexistence that 
disappear for = 0.40. The density and diffusion anomalous regions are also shifted to 
lower temperature, keeping the reduction in temperature-chemical potential phase diagram. 
The disappearance of the liquid-liquid temperature also reflects in the absence of density 
and diffusion anomalous regions in the limit of large density of obstacles. Both effects are 
related to both the entropy increase due to the presence of the obstacles and the disruption 
of the bonds network. 
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FIG. 11. Diffusion coefficient versus density at fixed temperature. The gray lines are fit of diffusion 
obtained by simulation. In case (A) we present the results for po = 0.08, in which the gray lines 
are just fit of diffusion obtained by simulation , blue dashed and dot-dashed connect the minimum 
and maximum in diffusion respectively. In (B) we have the same plot of case (A) but for the case 
Po = 0.24 and on inset we plot a zoom of anomalous region. In (C) the results for po = 0.40. The 
temperatures studied were T* = 0.30 ... 1.00 with AT* = 0.05 
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